In the previous section we explored spectra and developed an intuition for solving the peak pattern puzzle that each spectrum presents. Every x-ray fluorescence spectrum in our data cube is a sum of spectral patterns due to the presence of a mixture of chemical elements present at the measurement spot. These spectra are further complicated by 1) the presence of a continuum baseline caused by elastic and inelastic scattering of the paper background, and 2) Poisson photon noise. In other words, spectra contain mixed signals.

In order to arrive at element or material distribution maps, we now need to mathematically unmix the one million spectra from our data cube. Technically speaking we have to factorize our data cube $X$ into components $H$ and weights $W$(concentrations).

{% raw %} $$ X \approx H \times W $$ {% endraw %}

Different mathematical techniques have been developed for decomposing spectra into it's components. Best suited to our needs is a technique called Non-negative Matrix Factorization (NMF). Unlike other algorithms the NMF algorithm ensures that both the computed component spectra and their concentrations can not be negative.

The 2xNMF algorithm

The NMF algorithm is implemented in the scikit learn package and is fairly straightforward to use. In our case however we need to take into account that our data cube does not fit into memory all at once. As a workaround I will unmix/simplify the spectral data in two steps, applying the NMF algorithm twice.

STEP 1: Gaussian peak fitting

The idea of step one is resolve the problem of overlapping peak regions using NMF. This is done by first fitting Gaussian shaped peak profiles the all hotmax pixel peaks in all hotmax spectra. In order to get an idea what this means let's first fit a single spectrum using the fit_spectrum() function.

from maxrf4u import fit_spectrum 
from maxrf4u.peakmaps import _add_hotlines_ticklabels 

y = hotmax_spectra[5] # pick 
fitted_list = fit_spectrum(y, 'RP-T-1898-A-3689.datastack')

fig, ax = plt.subplots(figsize=[9, 5])

y_cont = fitted_list[-1]

for y_fit in fitted_list[:-1]: 
    ax.fill_between(x_keVs, y_fit + y_cont, y_cont, where=(y_fit>0.001), alpha=0.7)

ax.plot(x_keVs, y, label='spectrum to be fitted') 
ax.plot(x_keVs, y_cont, color='r', label='continuum')
ax.set_xlim([-1, 24])
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax);


ax.set_ylabel('Intensity (#counts)')
ax.set_xlabel('Energy (keV)')
ax.legend()
ax.set_title('Non-negative Matrix Factorization spectrum fit');

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-6f98119a30af> in <module>
     11 ax.plot(x_keVs, y_cont, color='r', label='continuum')
     12 ax.set_xlim([-1, 24])
---> 13 _add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax);
     14 
     15 

/mnt/datadisk/Work/Projecten/DoRe/code/maxrf4u/maxrf4u/peakmaps.py in _add_hotlines_ticklabels(datastack_file, ax, clip_vline)
    159     '''
    160     # read hotlines data
--> 161     ds = DataStack(datastack_file)
    162     x_keVs = ds.read('maxrf_energies')
    163     peak_idxs = ds.read('hotmax_pixels')[:, 2]

NameError: name 'DataStack' is not defined

FUNCTIONS

Compute continuum baseline from sum spectrum.

Uses rolling ball filter to remove peaks from sum spectrum.

Returns: y_continuum

Computes fitted and clipped Gaussian peak shapes for all hotmax pixels.

Returns: y_gauss_list

Fit Gaussian peak shapes to single spectrum `y`.

Uses NMF with fixed Gaussian component vectors.

Returns: weighed_components_list

import numpy as np 
import scipy.interpolate as sip
import scipy.signal as ssg 
from maxrf4u import DataStack, HotmaxAtlas 
from dask.diagnostics import ProgressBar 
import matplotlib.pyplot as plt 
import skimage.exposure as ske 
import skimage.morphology as skm 
from scipy.ndimage import gaussian_filter1d

y_continuum = get_continuum('RP-T-1898-A-3689.datastack')

y_gauss_list = get_gaussians('RP-T-1898-A-3689.datastack')

fig, ax = plt.subplots(figsize=[9, 5])
ax.plot(x_keVs, y_sum)
ax.plot(x_keVs, y_continuum)
for y_gauss in y_gauss_list: 
    ax.fill_between(x_keVs, y_gauss)

Focusing on Cu and Zn...

H = np.array([y_gauss for y_gauss in y_gauss_list]).astype(np.float32)
# add continuum 

H = np.r_[H, y_continuum[None,:].astype(np.float32)]

H.shape

(25, 4096)

n_components, _ = H.shape

X = hotmax_spectra[10:16]

W, H, n_iter = skd.non_negative_factorization(X, H=H, n_components=n_components, update_H=False)

X_nmf = W @ H

X_nmf.shape

(6, 4096)

W

array([[0.00556512, 0.04956888, 0.06137924, 0.23061307],
       [0.00457366, 0.1719793 , 0.03075115, 0.15105668],
       [1.0590024 , 1.1755654 , 1.2338187 , 0.75086987],
       [1.0590024 , 1.1755654 , 1.2338187 , 0.75086987],
       [1.0590024 , 1.1755654 , 1.2338187 , 0.75086987],
       [0.01936781, 0.10051632, 0.1715808 , 1.3440237 ]], dtype=float32)

nrows, _ = X.shape 
fig, axs = plt.subplots(nrows=nrows, sharex=True, figsize=[9, 12])

for i, y_nmf in enumerate(X_nmf): 
    axs[i].plot(x_keVs, X[i])
    axs[i].fill_between(x_keVs, y_nmf, color='orange')
    
axs[0].set_xlim([-1, 24])

(-1.0, 24.0)

       

# Plotting 

def plot_cube_slices(datastack_file, ax=None, tail_clip=0.05, xlim=[-1, 24]): 
    
    # read data 
    ds = DataStack(datastack_file) 
    x_keVs = ds.read('maxrf_energies')
    y_max = ds.read('maxrf_maxspectrum') 
    
    if ax is None: 
        fig, ax = plt.subplots(figsize=[9, 4]) 
        
    y_gauss_list = get_cube_slices(datastack_file, tail_clip=tail_clip)[1] 
    
    for y_gauss in y_gauss_list: 
        
        ax.fill_between(x_keVs, y_gauss) 
        
    ax.plot(x_keVs, y_max, color='r', alpha=0.5, label='max spectrum')
    ax.fill_between(x_keVs, y_max, color='r', alpha=0.2)
        
    _add_hotlines_ticklabels(datastack_file, ax)
    
    ax.set_xlim(xlim)
    
    ax.set_xlabel('energy [keV]')
    ax.set_ylabel('intensity [#counts]') 
    
    ax.legend()
    
    plt.tight_layout()
    
    return ax
        
        



# peak maps 

def get_peakmaps(datastack_file, slices, norm=True, verbose=False, algorithm='simple'): 
    '''Integrate peak `slices`  into peak maps and keV maps. 
    
    Returns: peak_maps, keV_maps'''
    
    if algorithm is 'simple': 
        print('Simply integrating peak slices without overlap correction')
    
    ds = DataStack(datastack_file)
    x_keVs = ds.read('maxrf_energies')
    cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
    
    peak_maps = []
    keV_maps = []
    
    # skipping baseline correction for now 
    # can not be based on super slices 
    
    n_channels = len(x_keVs) 
    
    print('Please wait two minutes while computing peak maps... ')

    for i, [si, sj, sk] in enumerate(slices): 

        peak_slice = cube[:,:,si:sk+1].compute() 

        # average over slice layers 
        print(f'Computing peak_map {i}/{len(slices)-1}...', end='\r')
        d = peak_slice.shape[-1]
        peak_map = np.sum(peak_slice, axis=2) / d 

        if norm:  
            peak_map = peak_map / peak_map.max()

        keV_idx_map = si + np.argmax(peak_slice, axis=2)
        keV_map = x_keVs[keV_idx_map]

        peak_maps.append(peak_map) 
        keV_maps.append(keV_map)
        
    print(f'Ready computing {len(slices)} peak maps.')
    
    return peak_maps, keV_maps  



def multi_plot(*images, hot_pixel=None, titles=None, roi_list=None, axis_off=False, 
               sharex=True, sharey=True, vmin=None, vmax=None, cmap='viridis', 
               fontsize='medium', zoom_xyc=None, zoom_half_wh=[100, 100]): 
    '''Inspect multiple images simultaneously... 
    
    Fold along multiple rows if n > 4'''
    
    nrows_max = 4 
    n_img = len(images)
    
    nrows = (n_img // nrows_max) # completely filled rows 
    rest = n_img % nrows_max 
    if rest != 0:
        nrows = nrows + 1
    
    if n_img <= nrows_max: 
        ncols = n_img
    else: 
        ncols = nrows_max
        
    figsize = [9, 5 + 1.3 * (nrows -1)]
    
    fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=figsize, squeeze=False, sharex=sharex, sharey=sharey)
    
    for i, img in enumerate(images): 
    
        axs.flatten()[i].imshow(img, cmap=cmap, vmin=vmin, vmax=vmax)
        
        if hot_pixel is not None: 
            hot_yi, hot_xi, hot_zi = hot_pixel 
            axs.flatten()[i].scatter(hot_xi, hot_yi, color='r')
        
        if zoom_xyc is not None:
            xc, yc = zoom_xyc
            w_ha, h_ha = zoom_half_wh 

            axs.flatten()[i].set_xlim(xc - w_ha, xc + w_ha)
            axs.flatten()[i].set_ylim(yc + h_ha, yc - h_ha)
        
        if roi_list is not None: 
            add_roi_patches(axs.flatten()[i], roi_list)
            
    if titles is not None:
        for i, t in enumerate(titles): 
            axs.flatten()[i].set_title(t, fontsize=fontsize)
    
    # seems to have no effect: 
    for i in range(n_img, nrows * ncols): 
        axs.flatten()[i].axis('off')
    # therefore trying this: 
    if axis_off: 
        axs_flat = axs.flatten()
        for ax in axs_flat: 
            ax.set_axis_off()
        
    fig.subplots_adjust(hspace=0.1, wspace=0.1)    
            
    plt.tight_layout()
    
    return fig, axs

from maxrf4u import plot_patterns, get_patterns, plot_cube_slices

ptrn_list = get_patterns(['S', 'Ca', 'K', 'Cl', 'Fe', 'Mn', 'Cu', 'Zn', 'Pb', 'Ti', 'Ba'])

fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 5])

plot_patterns(ptrn_list, ax=ax)
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax, clip_vline=False) 

plot_cube_slices('RP-T-1898-A-3689.datastack', ax=ax1); 


line_max, = ax1.plot(x_keVs, y_maxspectrum, color='r') 

ax1.set_ylim([-2, 70])


twax1 = ax1.twinx()
line_sum, = twax1.plot(x_keVs, y_sum, color=[0.3, 1, 0.3])

ax1.legend([line_sum, line_max], ['sum spectrum', 'max spectrum'])

ax.set_title('Peak pattern puzzle summary')
plt.tight_layout();

In order to ultimately create the much appreciated element distribution maps, we first need to slice the spectral image cube at all peak locations.

This can be done using the function get_cube_slices(). This function computes Gaussian peak shapes and corresponding slice indexes for each hotmax pixel in the spectral image cube. The result can be inspected with plot_cube_slices() function.

from maxrf4u import get_cube_slices, plot_cube_slices, get_peakmaps

slices, y_gauss_list = get_cube_slices('RP-T-1898-A-3689.datastack')

ax = plot_cube_slices('RP-T-1898-A-3689.datastack');
ax.set_title('Gaussian peak profiles')
plt.tight_layout()

If life were simple we could now integrate the intensities for a given peak slice to get the corresponding element map, and be done with our work! This is for example the case for our well-known, clearly visible and nicely isolated iron $K_{\alpha}$ peak located at 6.40 keV colored blue in the plot above. Let's check the precise limits of energy band that we want to integrate.

ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies') 

si, sj, sk = slices[10] # These are the Fe_Ka slice indices (left, mid, right) index  

print(f'The Fe_Ka peak located at {x_keVs[sj]:.02f} keV ranges from {x_keVs[si]:.02f} keV to {x_keVs[sk]:.02f} keV.')

The Fe_Ka peak located at 6.40 keV ranges from 6.04 keV to 6.76 keV.

Integration into an element map is simply reading the peak slice from the data cube, summing along the spectral axis and plotting the result.

# read Fe_Ka slice 
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
FeKa_slice = cube[:,:,si:sk+1].compute() # takes a few seconds... 
# integrate
FeKa_map = FeKa_slice.sum(axis=2)
# and plot...

fig, ax = plt.subplots(figsize=[8, 8])
ax.imshow(FeKa_map, vmax=75); # need to clip intensity due iron speckles  
ax.set_title('Simple peak slice map for iron Fe_Ka');

For now we just ignore errors due overlapping bands and take a preliminary look at all peak maps by simply integrating all peak slices use the get_peakmaps() function.

peak_maps, keV_maps = get_peakmaps('RP-T-1898-A-3689.datastack', slices, algorithm='simple')

Simply integrating peak slices without overlap correction
Please wait two minutes while computing peak maps... 
Ready computing 24 peak maps.

Although it is clear that these peak maps are not perfect yet, let's take a look. Starting point for the analysis of the maps is the summary of the peak pattern puzzle from the previous section.

Here is the overview of all 24 peak maps.

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

titles = [f'[{i}] {x_keVs[peak_idx]:.02f} keV' for i, peak_idx in enumerate(peak_idxs)]

fig, axs = multi_plot(*peak_maps_histeq, titles=titles)
fig.suptitle('Peak maps by simple integration')
plt.tight_layout()

Mixed signals

Unfortunately for most other peaks then iron XRF life is not so simple. Reason is that individual spectra contain mixed signals that originate from different chemical elements present at the pixel location where the spectrum was taken.

There are a two complications that we need to take into account. Especially for MA-XRF scans of drawings on paper we need to deal with many spectra for locations in contain only light elements. As a consequence a significant portion of their signal is due to inelastic (Compton) and elastic scattering by hydrogen atoms in the paper. We need to correct for this Compton ridge baseline when integrating these peaks. Our second problem is the occurrence of overlapping peaks. Let's first solve the baseline estimation problem...

The technical challenge here is that the baseline which is due to to scattering of the paper substrate will vary according to the local paper thickness, but also due to the local presence of heavier elements such as lead. A layer of lead atoms on top of a paper absorbs the incoming x-rays that therefore will not reach the paper and hence not be scattered to the detector. Spectra that contain lead peaks typically have a lower scattering baseline.

This behavior is exemplified if we compare two different hotmax spectra

The slowly varying ridge with Poisson noise and element specific emission peaks can clearly be seen in both spectra.

gap_slice, y_cont = get_continuum(x_keVs, y_sum, slices, plot_this=True)

Now I need to inspect how well this continuum baseline fits different spectra.

hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')

n = 4

y_hot = hotmax_spectra[n] 
y_snip = bsl.smooth.snip(y_hot, 50, )[0]

gap_slice, y_cont = get_continuum(x_keVs, y_sum, slices)

gap_average = get_gap_average(y_hot, gap_slice)

y_cont_baseline = y_cont * gap_average 

fig, ax = plt.subplots(figsize=[9, 4])

hma.plot_spectrum(n, ax)
ax.plot(x_keVs, y_cont_baseline, color=[0.3, 1, 0.3])
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax)
ax.set_xlim([-1, 24])
ax.set_ylim([-1, 20]);
#ax.plot(x_keVs, np.log(20*mu_Pb), color='k')
ax.plot(x_keVs, 0.6*y_maxsnip, color='orange')

[<matplotlib.lines.Line2D at 0x7f82d4281210>]

import xraydb

x_eVs = 1000 * x_keVs
mu_Pb = xraydb.mu_elam('Pb', x_eVs)

import pybaselines as bsl

y_maxspectrum

array([0.00239436, 0.00272245, 0.00339439, ..., 0.30186892, 0.31268886,
       0.31824136], dtype=float32)

y_maxsnip = bsl.smooth.snip(y_maxspectrum, 30, decreasing=True)[0]

y_maxsnip_noiseline = y_maxsnip + 0.1*y_maxsnip.max()

is_gap = y_maxsnip_noiseline > y_maxspectrum

fig, ax = plt.subplots(figsize=[9, 4])
ax.plot(x_keVs, y_maxspectrum)
ax.plot(x_keVs, y_maxsnip_noiseline)

ax.fill_between(x_keVs, y_maxsnip_noiseline, where=is_gap, color='orange', alpha=0.5)
ax.fill_between(x_keVs, y_maxsnip_noiseline, where=(labeled==19), color='r', alpha=0.5)
ax.set_xlim([-1, 25])
ax.set_ylim([-2, 50]);
ax.scatter(x_keVs[hill_top_idx], y_maxsnip[hill_top_idx])

<matplotlib.collections.PathCollection at 0x7f82d16d0150>

Need to inspect the width of each gap band to select best one. Or instead test dask boolean indexing.

is_test = np.zeros_like(x_keVs).astype(bool)

hill_top_idx = np.argmax(y_maxsnip)

is_test[1000:1100] = True

test_slice = cube[:,:,is_test]

test_slice

import skimage.morphology as skm

sizes[19-1]

27

labeled = skm.label(is_gap)

sizes = []
for l in range(1, labeled.max()): 
    
    d = np.sum(labeled == l)
    sizes.append(d)

sizes

[75,
 21,
 3,
 5,
 164,
 48,
 5,
 16,
 10,
 92,
 18,
 10,
 51,
 17,
 13,
 104,
 150,
 153,
 27,
 337,
 63,
 42]

y_sumsnip = bsl.smooth.snip(y_sum, 30, decreasing=True)[0]

fig, ax = plt.subplots(figsize=[9, 4])
ax.plot(x_keVs, y_sum)
ax.plot(x_keVs, y_sumsnip)

[<matplotlib.lines.Line2D at 0x7f82d3fa7950>]

slices, y_gauss_list = get_cube_slices('RP-T-1898-A-3689.datastack')

import matplotlib.cm as cm

y_continuum = smooth_compton(y_sum)

# peak free region for continuum normalization  
# not sure how to automate this choice 

# 18 - 19 is near continuum top and fairly large region 
left_peak_idx = 18 
right_peak_idx = left_peak_idx + 1
gap_band_i = slices[left_peak_idx][2] + 1
gap_band_j = slices[right_peak_idx][0] -1  

continuum_gap_mask = np.zeros_like(x)
continuum_gap_mask[gap_band_i:gap_band_j+1] = 1

print(f'Gap  region width: {gap_band_j - gap_band_i}')

Gap  region width: 105

y_

x = x_keVs 
y_maxspectrum = ds.read('maxrf_maxspectrum')
#def plot_cube_slices(slices, x, y): 

mask_list = []

for [si, sj, sk] in slices: 
    
    mask = np.zeros_like(x) 
    
    mask[si:sk+1] = 1
    
    mask_list.append(mask)
    
sum_mask = np.array(mask_list).sum(axis=0)
    
fig, ax = plt.subplots(figsize=[9, 4])

for i, mask in enumerate(mask_list): 
    
    ax.fill_between(x, 100*mask, where=mask, alpha=0.2)
    sj = slices[i][1]
    #ax.axvline(x[sj], color='r')
    
ax.plot(x, sum_mask)
ax.plot(x, y_maxspectrum) 
ax.plot(x, y_compton, label='Compton')
ax.legend()
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax)

ax.fill_between(x, y_maxspectrum, where=continuum_gap_mask, color='r')

ax.set_xlim([-1, 25]);

In this case, the best

datastack_file = 'RP-T-1898-A-3689.datastack'
tail_clip = 0.05 
xlim = [-1, 24]

#def plot_cube_slices(datastack_file, ax=None, tail_clip=0.05, xlim=[-1, 24]): 
    
# read data 
ds = DataStack(datastack_file) 
x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum') 

slices, y_gauss_list = get_cube_slices(datastack_file, tail_clip=tail_clip) 

n_slices = len(slices)

fig, ax = plt.subplots(figsize=[9, 4])

# tab20x2 color map 
# (there should be an easier way to cycle colors)
tab20 = cm.tab20(np.arange(20))[:,0:3]
tab20_r = cm.tab20_r(np.arange(20))[:,0:3]
tab20x2 = np.r_[tab20, tab20[2:]**0.4]
colors = tab20x2[0:n_slices]

# color gaussian peaks  

for i, y_gauss in enumerate(y_gauss_list): 

    ax.fill_between(x_keVs, y_gauss, color=colors[i]) 

# color corresponding slice  
for i, [si, sj, sk] in enumerate(slices): 
    
    is_slice = np.zeros_like(x_keVs)
    is_slice[si:sk+1] = 1  
    
    y_max = y_gauss_list[i].max() * np.ones_like(x_keVs) 
    
    ax.fill_between(x_keVs, y_max, where=is_slice, color=colors[i], alpha=0.3)

ax.plot(x_keVs, y_max, color='r', alpha=0.5, label='max spectrum')
#ax.fill_between(x_keVs, y_max, color='r', alpha=0.2)

_add_hotlines_ticklabels(datastack_file, ax)

ax.set_xlim(xlim)

ax.set_xlabel('energy [keV]')
ax.set_ylabel('intensity [#counts]') 

ax.legend()

plt.tight_layout()
    
# return ax

plot_cube_slices('RP-T-1898-A-3689.datastack')

<AxesSubplot:xlabel='energy [keV]', ylabel='intensity [#counts]'>

def plot_cube_slices(datastack_file, ax=None, tail_clip=0.05, xlim=[-1, 24]): 
    
    # read data 
    ds = DataStack(datastack_file) 
    x_keVs = ds.read('maxrf_energies')
    y_max = ds.read('maxrf_maxspectrum') 
    
    if ax is None: 
        fig, ax = plt.subplots(figsize=[9, 4]) 
        
    y_gauss_list = get_cube_slices(datastack_file, tail_clip=tail_clip)[1] 
    
    for y_gauss in y_gauss_list: 
        
        ax.fill_between(x_keVs, y_gauss) 
        
    ax.plot(x_keVs, y_max, color='r', alpha=0.5, label='max spectrum')
    ax.fill_between(x_keVs, y_max, color='r', alpha=0.2)
        
    _add_hotlines_ticklabels(datastack_file, ax)
    
    ax.set_xlim(xlim)
    
    ax.set_xlabel('energy [keV]')
    ax.set_ylabel('intensity [#counts]') 
    
    ax.legend()
    
    plt.tight_layout()
    
    return ax

Explorations

Zinc in the paper background

While trying to extract a smooth scattering baseline from the sum spectrum, my eye was drawn to a peak at the Zinc $K_{\alpha}$ energy. Because we see no copper, it seems that there is a low concentration of zinc present at a large surface in the drawing. I do not understand the source for that? Is is a modern conservation material? In order to find out I need a good peak map. Let's see if we can make one without yet having implemented the baseline estimation and peak overlap corrections...

fig, ax, twax = plot_cube_slices('RP-T-1898-A-3689.datastack')
ax.plot(x_keVs, 10 * y_sum, color=[0.3, 1, 0.3], label='sum spectrum')
ax.legend();

ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies') 

si, sj, sk = slices[13] # Zn_Ka slice (left, mid, right) index 

# read Zn_Ka slice 
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
with ProgressBar(): 
    ZnKa_slice = cube[:,:,si:sk+1].compute() # takes a few seconds... 
# integrate
ZnKa_map = ZnKa_slice.sum(axis=2)
# and plot

[########################################] | 100% Completed |  3.3s

fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])
ax.imshow(ZnKa_map, vmax=150); # need to clip?  
ax1.imshow(imvis, extent=extent)

x, y = 880, 555 
x2, y2 = 881, 559

ax.scatter(x, y, marker='+', color='r')
ax1.scatter(x, y, marker='+', color='orange')

ax.scatter(x2, y2, marker='+', color='r')
ax1.scatter(x2, y2, marker='+', color='orange')

<matplotlib.collections.PathCollection at 0x7f35dcb7ded0>

spectrum_xy = cube[x, y].compute()
spectrum_xy2 = cube[x2, y2].compute()

i, j, k = slices[13] # Zn 
Zn_band = np.zeros_like(spectrum_xy)
Zn_band[i:k+1] = 1

fig, ax = plt.subplots(figsize=[9, 4])
ax.fill_between(x_keVs, spectrum_xy, y_compton, where=(spectrum_xy>y_compton), label='xy1')
ax.fill_between(x_keVs, spectrum_xy2, y_compton, where=(spectrum_xy2>y_compton), label='xy2')


#ax.fill_between(x_keVs, spectrum_xy - y_compton, y_compton - y_compton, where=(spectrum_xy>y_compton), label='xy1')
ax.plot(x_keVs, y_compton, color='r')
ax.fill_between(x_keVs, 20 * Zn_band, where=(Zn_band>0), color='g', alpha=0.3)
ax.vlines([x_keVs[j]], 0, 20, color='g')
ax.legend()
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax)
mos.XFluo('Zn', tube_keV=40).plot(ax=ax);
ax.set_title('Hedgehawk plot for two bright Zn Ka pixels')
plt.tight_layout()

mos.XFluo('Zn', tube_keV=40).plot(ax=ax);

Compton peak

How does the gradient from left to right for zinc correlate to the Compton peak? I now wonder if the varying baseline in the cube is simply all related to the absorption of lead?

si, sj, sk = slices[21] # Rh_Ka Compton 

# read Compton_Rh_Ka slice 
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
with ProgressBar(): 
    Compton_slice = cube[:,:,si:sk+1].compute() # takes a few seconds... 
# integrate
Compton_map = Compton_slice.sum(axis=2)

[########################################] | 100% Completed |  6.3s

fig, [ax, ax1, ax2] = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=[9, 5])

ax.imshow(Compton_map, vmin=650); # need to clip?  
ax.set_title('Compton')
ax.scatter(*Compton_max[::-1], color='r')

ax1.imshow(imvis, extent=extent)
ax1.set_title('vis')

ax2.imshow(ZnKa_map, vmin=100, vmax=150); # need to clip?  
ax2.set_title('Zinc Ka');

Compton_max = np.argwhere(Compton_map==Compton_map.max()).flatten()
Compton_max

array([547, 127])

The 2xNMF algorithm

FUNCTIONS

`get_continuum`[source]

`get_gaussians`[source]

`fit_spectrum`[source]

Mixed signals

Explorations

Zinc in the paper background

Compton peak

The 2xNMF algorithm

FUNCTIONS

get_continuum[source]

get_gaussians[source]

fit_spectrum[source]

Mixed signals

Explorations

Zinc in the paper background

Compton peak

`get_continuum`[source]

`get_gaussians`[source]

`fit_spectrum`[source]